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Abstract 

We  present  a method  for  matching  a surface  in  three  dimensions  to  a set  of  data  sampled 
from  the  surface  by  means  of  minimising  the  distances  from  the  data  points  to  the  closest 
point  on  the  surface.  This  method  of  association  is  affine  transformation  invariant  and 
as  such  is  very  useful  in  situations  where  the  coordinate  axes  are  essentially  arbitrary. 
Traditionally,  this  problem  has  been  solved  by  minimising  the  (2  norm  of  the  distances 
from  the  data  points  to  the  corresponding  points  in  the  surface,  while  the  use  of  other 
tp  norms  is  less  well  known.  We  present  a method  for  template  matching  in  the  (\  norm 
based  upon  a method  of  directional  constraints  developed  by  Watson  for  the  related 
problem  of  orthogonal  distance  regression.  An  algorithm  for  this  method  is  given  and 
numerical  results  show  its  effectiveness. 

1 Introduction 

Template  matching  is  used  in  a variety  of  applications  such  as  the  quality  assurance  of 
manufactured  artifacts  [1]  and  dental  metrology  [2].  Given  a fixed  template,  i.e. , curve 
or  surface,  and  a set  of  data  in  a different  frame  of  reference,  template  matching  involves 
finding  the  frame  transformation  which  maps  the  data  onto  the  template. 

A typical  strategy  for  finding  the  optimal  transformation  parameters  in  the  template 
matching  problem  is  to  minimize,  in  some  norm,  the  orthogonal  distances  between  the 
transformed  data  and  the  template.  In  this  case,  the  template  matching  problem  can  be 
viewed  as  a form  of  orthogonal  distance  regression  (ODR)  [3],  which  is  a technique  com- 
monly used  for  fitting  curves  and  surfaces  to  measured  data.  Therefore,  most  algorithms 
for  solving  the  template  matching  problem  are  extensions  of  algorithms  for  ODR.  Tem- 
plate matching  in  the  £2  norm  is  addressed  by  Turner  [3]  and  in  the  £x,  norm  by  Butler 
et  al.  [1]  as  well  as  by  Zwick  [7]  for  the  two  dimensional  case. 

In  this  paper,  we  are  specifically  concerned  with  the  following  problem. 

Given  a fixed  differentiable  parametric  surface  f(u,u)  and  a set  of  m data 
{xJ-Tj  G 3i3,  find  points  {f(tq,Uj)}£L  l5  a rotation  matrix  Rq,  and  a transla- 
tion vector  t0  such  that  the  £\  norm  of  the  residual  distances  {||i?e(xj  - to)  - 
f(ui,Vi)\\2}iLi  is  minimal. 

This  is  the  template  matching  problem  in  the  £\  norm,  and  although  not  as  widely  used 
as  the  £2  and  £x  counterparts,  it  does  nonetheless  have  an  important  role  to  play.  The 
importance  of  the  £\  norm  is  that,  generally  speaking,  any  outlying  data  are  effectively 
ignored  with  the  result  that  an  approximation  is  obtained  which  is  largely  independent 
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of  any  unreliable  data.  This  has  particular  importance  when  our  data  arises  as  a result 
of  some  measurement  process,  perhaps  involving  many  complicated  and  finely-tuned 
instruments.  For  such  a measurement  scenario,  any  change  in  the  assumed  measurement 
conditions  can  result  in  a datum  which  has  gross  error  relative  to  other  data.  Thus,  if  we 
choose  a measure  which  is  susceptible  to  outlying  data,  we  are  in  danger  of  obtaining  an 
unrepresentative  approximation.  This  situation  is  avoided  by  use  of  the  l\  norm  and  we 
therefore  advocate  its  use  both  here  and  in  any  situation  involving  measurement  data 
where  a representative  approximation  is  required. 

A feature  of  optimal  t\  solutions  is  the  likelihood  of  a small  number  of  the  data  having 
a residual  of  zero,  and  it  is  therefore  unclear  whether  the  elements  of  the  Jacobian  matrix 
of  partial  derivatives  are  well-defined  for  these  points.  As  a result,  use  of  the  usual 
Gauss-Newton  method  would  appear  to  be  handicapped  due  to  its  dependence  upon 
the  Jacobian  matrix  to  calculate  an  updated  transformation  estimate.  This  difficulty 
also  arises  in  the  conventional  ODR  fitting  problem  and  has  recently  been  considered  by 
Watson  [6].  His  solution  is  to  adopt  a method  of  fitting  subject  to  directional  constraints. 
By  setting  these  directional  constraints  to  be  orthogonal  to  the  approximant,  Watson 
shows  not  only  that  the  Jacobian  is  defined  but  also  how  to  compute  its  elements  without 
incurring  a build-up  of  rounding  error. 

In  this  paper,  we  extend  Watson’s  constrained  direction  fitting  routine  to  the  template 
matching  problem.  We  show  that  Watson’s  results  are  equally  valid  for  iy  template 
matching.  Finally,  we  exploit  these  results  to  give  a reliable  algorithm  for  the  iy  template 
matching  problem. 

The  structure  of  this  paper  is  as  follows.  Section  2 provides  the  results  necessary 
to  justify  the  new  technique.  Section  3 describes  the  algorithm  adopted  to  implement 
the  theory.  Section  4 gives  some  numerical  results  for  both  a simple  case  and  a more 
challenging  case.  Finally,  Section  5 concludes  this  paper  and  presents  possibilities  for 
future  work. 

2 Theory 

We  are  concerned  with  the  minimisation  of  the  quantity 

E = ||(d1,...,dm)||,  (2.1) 

where 

di  = min  ||x;  — i(ui,vi)\\2  , i = 1,2,.'..  ,m,  (2.2) 

U{,Vi 

and 

x = f?e(x- 10),  (2.3) 

with  respect  to  the  rotation  parameters 
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the  translation  parameters 

to 

and  the  location  parameters 

U 

This  is  a constrained  problem  and  can  be  solved  using  a separation-of-variables  approach 
as  described  by  Turner  [3]  among  others.  In  this  approach,  the  problem  of  obtaining  the 
transformation  parameters 

-(£)■ 

is  separated  from  the  subproblem  of  obtaining  the  location  parameters  U.  At  each  iter- 
ation, the  subproblem  is  solved  to  obtain  an  optimal  U for  the  current  transformation 
parameters  t which  is  then  used  to  obtain  an  update  of  the  transformation  parameters 
themselves. 

2.1  Considerations  specific  to  the  t\  problem 

Up  to  this  point,  we  have  not  specified  which  norm  we  are  using  to  measure  the  disparity 
between  the  transformed  data  and  the  template.  Since  we  will  be  particularly  interested 
in  the  i\  case,  this  section  discusses  problems  inherent  in  the  solution  of  such  a problem. 

The  major  problem  with  solving  non-linear  problems  is  that  in  order  to  use  a 
technique  such  as  the  Gauss-Newton  method,  derivative  information  is  required.  Unfor- 
tunately, derivatives  of  the  distances  d are  not  defined  when  a distance  has  a value  of 
zero.  Such  is  the  nature  of  i\  approximation  that  zeros  are  to  be  expected  at  an  optimal 
solution  [5],  Thus,  it  is  unclear  whether  the  Jacobian  matrix  is  defined  at  these  data 
points.  Recent  work  by  Watson  [6]  has  considered  how  the  related  problem  of  ortho- 
gonal distance  regression  might  be  solved  by  considering  distances  to  be  measured  along 
fixed  direction  vectors  w;.  Orthogonal  distance  regression  involves  the  fitting  of  a curve 
or  surface  to  a set  of  data  where  the  residuals  are  taken  to  be  the  shortest  distance  from 
the  data  to  the  approximant  [3].  Template  matching  can  be  seen  as  a variant  of  this  since 
the  residuals  are  measured  in  the  same  way,  but  we  are  only  altering  the  position  and 
orientation  of  the  approximant,  rather  than  the  actual  shape  itself.  Thus,  techniques  for 
orthogonal  distance  regression  can  be  used  successfully  in  template  matching. 

By  means  of  these  directional  constraints,  it  is  possible  to  show  that  if  we  choose  the 
directions  w*  to  be  the  orthogonal  directions, 

f(Wj,Uj)  Xi 

w — i i 

||f(w*,«i) -X»l|2 

then  the  derivatives  are  well  defined  in  the  limit  as  ||f(uj,u,)  - 5c, 1 1 2 — » 0. 

This  result  may  be  summarised  in  the  following  Theorem  (taken  from  Watson  [6]). 
Theorem  2.1  For  parametric  fitting,  let  the  (usual)  Gauss-Newton  method  produce  a 
sequence  {t}  such  there  is  a unique  unit  normal  vector  to  the  template  at  f («,■ , vt ) , and 


Template  matching  in  the  i\  norm 


133 


x,  remains  on  one  side  of  the  template.  Then  Vt di  is  well  defined  on  this  sequence. 

If  f (ui,Vi)  — > x,,  then  this  formulation  will  lead  to  similar  problems  to  which  we 
are  attempting  to  resolve  as  a result  of  the  quotient  becoming  undefined.  As  a result, 
Watson  [6]  suggests  leaving  w*  unchanged  once  di  becomes  small.  By  this  method, 
numerical  problems  arising  as  a result  of  a distance  tending  to  zero  may  be  avoided. 
However,  the  algorithm  will  still  tend  to  the  correct  solution  provided  that  the  small 
residual  corresponds  to  an  interpolation  point  of  the  l\  solution.  If  this  is  not  the  case, 
then  the  solution  will  not  be  optimal,  but  will  still  be  close  to  the  optimal  solution. 

2.2  Possible  problems 

The  most  immediate  problem  that  arises  is  how  to  ensure  that  there  exists  a point  on 
the  template  which  is  situated  along  the  direction  vector  given  from  each  datum.  Clearly 
in  certain  situations,  there  will  not  exist  such  a point  — corresponding  to  the  case  where 
the  direction  vector  lies  within  the  tangent  plane  of  the  template  in  the  region  of  the 
datum.  In  such  a situation  there  would  seem  to  be  two  possible  recourses  available. 

(1)  Ignore  these  data. 

(2)  Choose  the  point  on  the  template  that  is  closest  to  the  line  though  the  datum 
defined  by  the  direction  vector. 

It  has  been  found  through  empirical  results  that  provided  the  problem  only  occurs  on 
certain  iterations  rather  than  as  a result  of  poor  choice  of  the  direction  vectors  associated 
with  the  template,  ignoring  the  problem  data  is  the  better  option.  Use  of  the  second 
option  has  been  found  to  prevent  convergence  of  the  algorithm. 

3 Algorithm 

The  algorithm  to  implement  this  technique  consists  of  two  sub-algorithms,  each  related 
to  a specific  section  of  the  main  algorithm.  These  sub- algorithms  are 

(1)  the  constrained  closest  point  problem, 

(2)  the  calculation  of  a new  transformation  estimate. 

3.1  Constrained  closest  point  problem 

For  each  data  point  x;,  this  problem  is  that  of  finding  Ui  and  v.L  such  that  the  constraint 

x — i(u,  v)  = dvr,  (3.1) 

is  satisfied  (subscripts  dropped  for  clarity).  Expanding  this  equation,  we  obtain 


If  we  pre-multiply  this  equation  by  aT,  we  obtain 

aTx  — aTf(u,  v)  — daTw  = 0.  (3.2) 
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Thus,  by  choosing  a to  be  orthogonal  to  w,  we  are  able  to  eliminate  d from  equation 

(3.2) .  Similarly,  if  we  multiply  equation  (3.1)  by  b we  obtain  the  equation 

b'  x - bTf  («,?.')  - dbTw  = 0. 

We  may  thereby  reduce  the  system  (3.1)  to  that  of  two  (nonlinear)  equations  in  two 
unknowns  ( u and  v).  This  system  can  then  be  solved  by  adopting  a Newton- type  method. 
Our  problem  has  been  reduced  to  that  of  solving 

F(u,v)  = [a  : b]T(x  - f(u,v))  = 0, 

which  has  derivative 

VU)„F  = —[a  : b]T(V„f  : Vt,f), 

by  means  of  Newton’s  method  which  involves  adopting  an  iterative  approach  and  solving 

V"'”F(  SSv)  = ~FM  (3’3) 

at  each  stage  to  obtain  better  estimates  u + 5u  and  v + Sv.  The  quantities  F(u,  v)  and 
VUtVF  are  straightforward  to  calculate  as  they  arise  directly  from  the  explicit  paramet- 
erisation  of  the  template. 

All  that  remains  is  the  choice  of  a and  b.  We  obtain  these  vectors  by  taking  the  cross 
product  of  w with  two  arbitrary  vectors  — resulting  in  two  vectors  which  are  orthogonal 
to  v.  More  generally,  the  vectors  a and  b should  be  chosen  to  ensure  that  the  system 

(3.3)  is  well-conditioned. 

3.2  Updating  the  transformation  estimate 

The  method  we  adopt  to  obtain  an  update  of  the  transformation  parameters  is  the 
Gauss-Newton  method.  This  involves  solving,  at  each  iteration,  the  problem 

«7<5t  = — d,  (3.4) 

in  the  l\  sense,  where  J is  the  Jacobian  matrix  of  partial  derivatives  with  entries  Jy  = 
Vt3.dj.  The  estimate  of  the  optimal  transformation  parameters  is  then  updated  according 
to 

t = t + <5t. 

Thus,  since  the  distances  d are  obtained  from  the  constrained  closest  point  subprob- 
lem, we  are  left  with  the  task  of  calculating  the  Jacobian  matrix.  For  each  datum,  from 
equation  (3.1),  we  have  that 

x(t)  - f(u(t),u(t))  = wd(u(t),  v(t)), 

where  we  have  explicitly  included  the  dependency  of  the  distance  d on  the  location 
parameters  U.  Differentiating  and  rearranging,  we  obtain 

Vtx  = wVtd  + V[/fVtf/. 

This  is  equivalent  to  the  form 

Vtx  = [w:V,;f](  ^ 


Template  matching  in  the  l\  norm 


Therefore, 

J = Vtd  = ef[w  : Vuf]_1Vtx, 

where  ei  is  the  first  component  vector.  Having  obtained  the  Jacobian  matrix  J and  the 
distance  vector  d,  we  are  now  in  a position  to  solve  the  system  (3.4)  in  order  to  update 
our  estimate  of  the  optimal  transformation  parameters  t. 

We  note  that  using  the  traditional  orthogonal  distances  can  lead  to  problems  since 
calculation  of  the  Jacobian  matrix  involves  division  of  each  row  by  the  corresponding 
orthogonal  distance  — leading  to  exacerbation  of  rounding  errors  and  possible  division 
by  zero  especially  in  the  case. 

4 Numerical  results 

In  this  section,  we  present  two  example  to  illustrate  the  techniques  presented  in  this 
paper.  In  the  first,  we  have  a small  number  of  data  which  we  wish  to  match  to  a given 
plane.  In  the  second,  we  have  a larger  number  of  data  and  we  wish  to  match  them  to 
a cylinder.  In  both  cases,  although  analytical  expressions  are  available  to  obtain  the 
constrained  closest  points  on  the  templates,  we  nonetheless  utilise  the  method  presented 
above  in  order  to  test  its  effectiveness. 


4.1  Simple  problem 


Here  we  describe  the  problem  of  matching  a representative  set  of  8 data  onto  the  plane 
defined  as 


f (it,  v)  = u 


1 

0 

0 


+ v 


0 

1 

0 


Since  this  problem  is  rank  deficient  if  we  use  all  six  possible  transformation  parameters, 
we  restrict  ourselves  to  using  a translation  in  the  ^-direction  and  rotations  about  the  x 
and  y axes. 

Having  three  degrees  of  freedom,  we  might  expect  to  obtain  an  optimal  ti  solution 
which  interpolates  3 of  the  data.  However,  as  we  shall  see,  this  is  unattainable  in  general 
and  we  can,  in  fact,  only  expect  interpolation  at  two  points.  As  Watson  states  [6],  in 
such  a situation,  the  rate  of  convergence  can  be  unacceptably  slow.  This  is  found  to  be 
the  case.  It  can  be  seen  that  not  only  is  the  convergence  slow,  but  an  optimal  solution 


Iteration 

norm  (residuals) 

norm  (update) 

1 

0.6662 

4.9901e-02 

5 

0.3008 

3.5716e-04 

10 

0.3007 

8.8545e-06 

50 

0.3006 

9.1533e-06 

100 

0.3008 

3.8514e-04 

Tab.  1 Progress  of  the  Gauss-Newton  method  for  planar  data. 


is  never  obtained,  with  the  objective  function  ||d||i  increasing  occasionally. 
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To  ensure  convergence,  a simple  line-search  algorithm  was  adopted  which  searches 
along  the  direction  obtained  from  the  Gauss-Newton  step  for  the  maximum  reduction 
in  the  objective  function.  This  modification  affects  convergence  in  3 iterations. 

4.2  A more  challenging  problem 

As  a more  challenging  problem,  we  consider  the  matching  of  a set  of  128  data  which 
supposedly  represent  a cylinder  but  which  contain  8 wild  points.  The  cylinder  is  para- 
metrised by  u and  v as 

(cosu  \ 
sinu  , 

v I 

resulting  in  a cylinder  with  unit  radius  oriented  along  the  z-axis.  Again,  the  problem 
of  matching  the  data  onto  this  model  is  rank  deficient.  The  rank  deficiencies  occur  due 
to  rotations  about  the  z-axis  and  translations  along  the  z-axis.  As  such,  we  omit  these 
possible  transformations. 

Although  we  might  initially  expect  to  interpolate  4 data  points  at  an  optimal  solu- 
tion, we  find  that  in  fact  only  two  are  guaranteed,  although  if  a third  point  lies  within 
two  radii  of  one  of  these  two  points,  then  three  points  can  be  guaranteed.  Typically,  this 
will  occur  when  the  data  is  representative.  For  the  data  set  we  are  considering,  we  expect 
three  interpolation  points  due  to  the  data  representing  the  cylinder  and  in  fact  at  the 
optimal  solution,  three  interpolation  points  are  obtained.  In  fact,  the  “missing”  interpol- 
ation has  the  effect  of  slowing  convergence  of  the  Gauss-Newton  method  considerably  so 
that  in  100  iterations,  the  algorithm  had  not  been  deemed  to  converge.  However,  by  the 
introduction  of  a simple  line-search  method,  the  algorithm  converged  in  five  iterations 
as  displayed  in  Table  2. 


Iteration 

norm(residuals) 

norm(update) 

1 

0.9654 

5.6796e-03 

2 

0.9559 

6.2932e-04 

3 

0.9557 

1.0141e-04 

4 

0.9557 

2.5812e-07 

5 

0.9557 

4.4006e-14 

Tab.  2 Progress  of  the  Gauss-Newton  method  for  cylindrical  data  using  a line-search. 


5 Conclusions 

This  paper  has  shown  how  perceived  problems  in  (i  template  matching  can  be  avoided 
by  use  of  the  so-called  “method  of  directional  constraints”.  In  this  method,  the  closest 
point  on  the  template  along  a given  direction  vector  is  calculated  in  order  to  obtain  the 
residuals  between  data  and  template.  By  then  altering  this  direction  vector  to  be  the 
normal  to  the  surface  at  that  projected  point,  the  algorithm  progresses  to  the  expected 
i\  solution.  Problems  regarding  undefined  quotients  are  avoided  by  no  longer  updating 
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the  direction  vectors  corresponding  to  a datum  when  the  residual  associated  with  that 
point  is  below  a certain  tolerance. 

This  work  forms  part  of  a larger  project  to  consider  novel  approaches  to  ill-conditioned 
problems  in  metrology.  It  is  hoped  that  the  work  presented  in  this  paper  will  aid  in  the 
resolution  of  rank-deficient  systems  and  ill-conditioned  systems  by  altering  the  usual 
orthogonal  distances  to  be  these  directional  constraints,  which  should  remove  some  of 
the  rank  deficiency. 

As  an  example,  consider  the  template  matching  problem  where  the  template  to  be 
matched  is  an  infinite  cylinder  with  axis  along  the  2-axis.  Using  typical  template  match- 
ing algorithms,  this  problem  is  rank  deficient  by  two  at  the  solution  due  to  the  possible 
translation  in  the  2-axis  and  the  possible  rotation  about  the  2-axis.  By  introducing  these 
directional  constraints,  the  rotational  rank  deficiency  is  almost  completely  resolved  (there 
are  now  two  possible  rotations  to  obtain  the  optimal  matching  rather  than  the  infinite 
number  previously). 

The  use  of  the  i\  norm  is  also  being  used  to  attempt  and  resolve  any  rank  deficien- 
cies and  ill-conditioning  present  in  the  problem.  This  is  achieved  by  ensuring  that  any 
local  deviations  from  the  template  (caused  by,  for  example,  wear)  are  “ignored”  so  that 
regions  of  local  deviations  might  be  compared.  This  will  then  result  in  a resolution  of 
the  uncertainty  in  the  transformation  parameters. 
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